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abstract 

When  electron-electron  correlations  are  important,  it  is  often  necessary  to  use  "e.vact”  numerical 
methods,  such  as  Lanczds  diagonalization,  to  study  the  full  many-body  Hamiltonian.  Unfortunately, 
such  exact  diagonalization  methods  are  restricted  to  small  system  sizes.  We  show  that  if  the  Hubbard 
U  term  is  replaced  by  a  “periodic  Hubbard”  term,  the  full  many  body  Hamiltonian  may  be  exactly 
solved,  even  for  very  large  systems,  though  for  low  fillings.  However,  for  half-filled  systems  and  large 
U  this  approach  is  not  only  no  longer  exact,  it  no  longer  improves  extrapolation  to  larger  systems. 
We  discuss  how  generalized  “randomized  variable  averaging”  (RVA)  or  “phase  randomization”  schemes 
can  be  reliably  employed  to  improve  extrapolation  to  large  system  sizes  in  this  regime.  This  general 
approach  can  be  combined  with  any  many-body  method  and  is  thus  of  broad  interest  and  applicability. 

INTRODUCTION 

Solutions  of  the  Peierls- Hubbard  Hamiltonian  (PHH)  (1)  for  small  clusters  are  strongly  dependent 
on  the  boundary  condition  (BC);  z.c.,  for  the  case  of  polyacetylene,  whether  the  carbon  atoms  are 
viewed  as  being  on  a  real  ring,  such  as  benzene,  a  real  chain,  such  as  (l,3,5)-hexatriene,  or  a  more 
exotic  geometry,  such  as  anti-periodic  BCs.  We  have  found  that  (2,3)  to  reduce  finite  size  corrections 
to  the  calculated  optical  gaps  and  spectra,  thus  improving  the  extrapolation  to  the  infinite  case,  it  is 
effective  to  employ  novel  “averaging”  techniques  to  “randomize”  the  many-body  levels  which  influence 
those  correlation  functions  of  interest.  Briefly,  by  randomization,  we  mean  simply  that  the  total  energy 
of  the  system  is  to  be  viewed  as  a  weighted  average  of  the  energies  derived  separately  for  each  of  several 
different  “random”  values  (RV)  of  a  particular  parameter,  which  may  be  the  BC,  a  hopping  integral, 
an  on-site  energy,  or  the  Hubbard  U:  E  =  where  xrv  is  a  normalized  weighting 

factor  with  =  !•  Similarly,  any  correlation  function  can  be  viewed  as  the  average  of  of 

its  values  from  the  separate  RVs.  The  number  of  RV,  Wrv,  may  be  varied,  and  the  choice  of  the  RV 
and  zRv  may  be  truly  random,  or  may  be  a  pre-specified  set  of  values,  which  may  in  turn  depend  on 
the  quantity  being  studied.  The  case  of  phase  randomization  and  its  rdation  to  Bloch's  theorem  are 
described  bekm.  More  ge’^rrally,  this  randomized-variable-averaging  (RVA)  technique  may  be  viewed 


t 


We  lix  us  in  this  manuscript  on  tlio  calculation  via  exact  Lanczos  diagonmization  of  tlic  self- 
consistent  uniform  dimerization  and  one-photon  optical  absorption  of  the  ground  state  within  the  l-I). 
1/2-filled,  one-band  PHH.  However,  RVA  is  equally  applicable  to  the  calculation  of.  e.g.,  the  phonon 
modes  of  the  ground  state,  as  well  as  the  geometry  and  self-consistent  absorptions  of  doped  or  higher 
lying  (e.^.,  triplet)  states  [3],  luminescence  spectra,  structure  factors,  etc.,  in  single-  or  multi-band 
models  in  1-,  2-,  and  3-D,  via  Lanczos  diagonalization,  as  well  as  other  numerical  procedures  on  small 
lattices,  such  as  Monte  Carlo. 


PHASE  BOUNDARY  CONDITION  AVERAGING  AND  BLOCH’S  THEOREM 

Bloch’s  theorem  tells  us  that  the  single  particle  wavefunctions  of  a  system  of  size  M-N  with 
periodicity  N  and  periodic  boundary  conditions  (PBCs)  are  exactly  the  single  particle  wavefunctions 
of  the  ensemble  of  systems  of  size  N  for  each  of  the  phase  BCs  -f  1)  =  e*‘’c'’(l),  4)  =  2irifM. 
I  =  1,...,M.  Restated,  the  (exact)  properties  of  a  larger  (single  particle)  system  may  be  found  by 
forming  a  (particular)  average  over  smaller  systems  with  different  BCs.  This  "boundary  condition 
averaging”  (BCA)  or  “phase  randomization”  is  a  special  case  of  our  general  RVA  procedure.  Blochs 
theorem  may  be  generalized  to  many  particle  wavefunctions  and  is  applicable  here  if  the  form  of  the 
Hamiltonian  is  modified  slightly.  The  one-dimensional  PHH  is 

H  =  ^(-<0 +a<5/)(cj^c/+i  ^ -f  <,)  -I-  (1) 

i,<r  t  ^  t 

If  we  replace  the  Hubbard  term,  U  'ts  periodic  analog, 

U  ^  ^ 

Hel-el  ^  ,  (2) 

”  tsl 

then,  if  the  original  problem  had  periodicity  N ,  {ne,a)  =  {nt+vN.a),  we  have,  in  a  mean-field  sense, 
the  “same”  Hamiltonian.  With  this  new  “periodic-Hubbard”  Hamiltonian  (periodic  PHH)  the  Bloch 
analysis  of  a  given  large  (periodic)  system  {N  -M  sites,  Ng-M  electrons)  is  accomplished  by  considering 
a  small  system  (N  sites,  Ng  electrons)  with  several  different  BCs.  The  method  scales  linearly  with  M, 
allowing  one  to  handle  reasonably  large  systems,  and  sets  a  sound  theoretical  basis  for  the  empirical 
observation  [2]  that  RVA  can  be  used  to  smooth  optical  absorption  spectra  obtained  via  exact  finite-size 
diagonalization. 

Using  the  periodic-Hubbard  term,  Eq.  (2),  assuming  the  lattice  distortion  has  the  same  periodicity, 
6{£  +  N)  =  6{£),  and  using  2<7=i.i  =  Ng-  M,  one  can  show: 


M  S 


tt.vssl  i,j,k,l=l 


where 


ff(i,j,k,l;si,S2)  =H{i  -f-  N,j  +  N,kJ;si,S2) 

—  ^(f,  J,  1,  ^2,  Si  )  —  ff(fc, /,  t,  J  ,  Sj ,  S2  )  —  H k,  1]  S\,  S2) 


1  r  ^  a 

H{i,j,k,l;si,S2)  =  -  ^(-fo  +  a5„)(3i,„5j,„+i -h  +  {ij  ^  k,l) 

^  '•„=i 

n=l  *  * 


Note  that  H(i,j^k,l;s\,S2)  is  independent  of  fi,  and  A/.  I'he  Uloch  analysis  on  the  many-body 
eigenfunctions  using  the  symmetries  of  H{i^j,k,l-,Si,S2)  leds  to  a  of  the  form; 


cr,,...,cr,  =  t,i  ni,...,n,=l 


+  lipN,a. 


where  V'fc, . kt . n<,  ...,np)  =  V»*, . . *,(«!,  -,  ni+fn-N, tip) ,  ke  =  (2irj/M),  je{0, M 

1},  and  p  =  Ne  •  M .  We  find  for 


N  N 


ki,...,kp  —  ^  ^  ^  ]  ^k, . 1  ^p)  ,j  C  * 


<Ti,...,«r,=  {,i  ni . np  =  l  e=\ 


AWM  tWB 


(i/o  comes  into  play  for  operators  like  c)y,  ,CAr+i,j).  From  Eq.  (5)  we  see  Rper  is  diagonal  in  the  It’s, 


. k,  =  . '  ft)! 

as  are  p  =  Yl,i,a  ^t,<r  J  =  *  12t,c(~^o  +  ot^e)ict^„ce+i^„  -  and  the  eigenfunctions  of  Hper 

can  be  written  in  the  form  Eq.  (4). 

The  symmetry  of  9  implies  that  if  a/,  =  at,  and  =  kt,  then  nr,  ^  nt, .  Thus  there  are  at  most 
2N  of  the  kt  the  same.  We  postulate  that  the  ground  state  lies  in  the  manifold  with  each  of  the  M 
distinct  values  for  kt  occurring  Ng  times  and  (for  the  half- filled  bzmd)  with  equal  numbers  of  up  and 
down  spins  (in  general,  the  largest  manifold).  This  can  be  checked,  and  we  stress  that  the  analysis  up 
to  this  point  is  exact  for  the  periodic  PHH.  In  this  manifold  we  can  write  as: 

2N 

=  Y  •••  ,  (7) 


aj.— .or»f  =  1 


where  qt  =  {2Trl/M)  and 


=  E  E  E 

<^1  =  T  . «'/»,  =  !  »»i.— 

M  M 

X  V  ...  V  /.f  mn, 

<*1  =  1  <*w,  =1 

where  N\  =  7V|  =  [[^]]  ([( ]]  denotes  int^er  value  of).  Note 

ci+yv.<,civ.<^  =  «”  Ci.<,CAr,„  , 


and  so  solving  the  Af*  electron  problem  on  the  full  N  ‘M  sites  with  PBCs  for  fixed  q  is  exactly  equivalent 
to  solving 

^per  (9) 

on  N  sites  with  the  q-dependent  BC  defined  by  EJq.  (8).  Thus,  to  find  the  exact  Ng  electron  eigenfunc¬ 
tions  on  N  •  M  sites,  we  only  need  to  solve  Eq.  (9)  with  N  sites  and  N,  electrons  for  each  of  the 
M  BCs,  rather  tham  the  full  N  -  M  site  problem. 


To  obtain  the  exact  eigenfunctions  of  the  full  jV^  A/  electron  problem,  u.'  need  to  .solve  for 
(//’(oi, ...,  qm)-  This  is  where  an  approximation  must  be  made  to  be  able  to  solve  the  problem  numer¬ 
ically,  as  we  assume  N  is  already  as  large  as  computationaUy  feasible.  The  full  eigenvalue  problem 
yields 

N  ,  r  M  -1  Af 


N. 

<l.a 


= 


N  f  r  M  Af 

^  /=1  -I  i/=I 


V.l 


f  M 

5 


nt+uN,i  1  ^ 


N. 

9.0 


M  . 


(10) 


To  solve  Eq.  (6)  we  must  use  some  approximation  for  Eq.  (10),  such  as  perturbation  theory  or  mean- 
field.  Thus  BCA  means  we  treat  electron-electron  correlations  within  the  Ng  electron  manifold  exactly, 
and  electron-electron  correlations  between  Ng  electron  manifolds  approximately.  The  BCA  results 
reported  here  are  zero-order  perturbation  theory:  we  have  assumed  ^(oj , ...,  qay )  to  be  a  product  of 
A-functions;  i.e.,  for  the  ground  state  we  use: 


M 


^91.0 


*9«.0  ’ 


t=l 


(11) 


Note  that  if  we  had  been  interested  in  systems  at  low  filling  (Ng  electrons),  rather  than  near  half  filling 
(Ng  •  M  electrons),  there  would  be  no  need  to  make  any  approximations. 


Comparison  of  the  periodic  and  standard  Peierls-Hubbard  Models  at  half-filling 

At  U =0  Bloch’s  theorem  for  the  many  particle  wavefunctions  as  formulated  above  is  exau:t  at  all 
fillings.  For  the  half-filled  PHH,  we  have  just  shown  that  phase  BCA  involves  two  approximations: 
first,  replacing  the  Hubbard  term  by  Eq.  (2),  and,  second,  ignoring  correlations  between  the  Ng  particle 
wavefunction  leading  to  Eq.  (11).  We  now  test  these  approximations.  In  Fig.  1  we  show  the  minimum 
energy  dimerization  of  an  8-site  system  as  a  function  of  U  for  (t)  phase  averaged  solution  and  (tt) 
the  exact  solution  of  the  half-filled  band  with  the  periodic  Hubbard  term  and  periodicity  2,4,  and 
8  (periodicity  8  is  the  the  usual  Hubbard).  The  approximate  periodic  (BCA)  solution  lies  between 
the  exact  periodic  solution  and  the  exact  pure-Hubbard  solution.  Since  the  originad  problem  was  for 
U  •  Af  53^1  on  N  ■  M  sites,  one  might  argue  that  the  approximate  phase-averaged  result  for 

^  SrLi  Af  sites  may  be  closer  to  the  desired  answer  than  the  exact  periodic-Hubbard  N  •  M 

site  result. 

In  Fig.  2  we  show  the  dimerization  amplitude  obtained  by  phase  BCA.  We  see  that  the  extrapolated 
infinite  behavior  at  small  to  intermediate  U  is  well  approximated  after  phase  averaging  for  N=10  even 
without  extrapolation.  The  agreement  is  best  for  small  U,  while  no  change  from  the  pure  periodic 
behavior  is  found  for  large  U.  That  this  particular  technique  is  expected  to  have  no  effect  at  large  U 
cam  be  seen  by  examining  the  effective  spin-Peierls  Hamiltonian  [4],  which  is  asymptotically  independent 
of  the  phase  of  the  BC.  However,  the  phase-averaged  analysis  does  lead  to  the  interpretation  of  the 
finite-size  result  as  the  result  for  the  larger  system  with  an  effective  periodic-Hubbard  interaction. 
Thus,  if  U  is  scaled  by  the  system  size,  as  is  done  in  Fig.  2,  the  infinite  behavior  should  be  more  easily 
extrapolated.  It  appears  that  this  is  at  least  approximately  true. 

To  show  how  phase  BCA  affects  the  optical  spectrum,  in  Fig.  3a  we  show  the  (Lorentzian  broad¬ 
ened)  spectrum  at  ff/to=0.4  for  N=8  and  5  phase  BCs  (enough  to  yield  a  smooth  spectrum  at  U=0). 
It  is  clear  how  the  different  BCs  “fill  in”  the  spectrum.  Fig.  3b  shows  the  spectrum  nt  U/to=l.6.  Here, 
despite  the  several  BCs  used,  the  spectrum  remains  sp2irse,  due  to  the  loss  of  effectiveness  of  the  phase 
BCA  scheme  as  discussed  above.  If  we  were  to  use  a  broad  enough  Lorentzian  to  smooth  this  spectrum 
(the  same  width  as  in  Fig.  3a  was  used),  the  optical  gap  edge  would  be  lost. 
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Fig.  1  (left).  Comparison  of  the  approximate  and  exact  periodic-Hubbard  minimum  energy  dimerization 
as  a  function  of  U  for  to=2.5  eV,  a=4.1  eV fk,  and  K=2l  eV/A^. 

Fig.  2  (right).  Dimerization  vs.  U  for  various  N:  phase  BCA.  Parameters  as  Fig.  1. 
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Fig.  3.  The  Lorentzian  broadened  optical  absorption  spectra  calculated  using  the  phase  BCA  method 
for  (a)  U=1  eV  and  (b)  U=4  eV.  Note  the  “sparseness”  of  the  spectrum  in  the  latter  case  indicating 
the  failure  of  the  phase  BCA  approach.  For  comparison,  we  also  show  the  spectrum  at  [7=10  eV  using 
(c)  amplitude  BCA  and  (d)  an  RVA  procedure  where  U  is  varied  on  one  site.  Note  these  RVA  spectra 
are  dense  and  show  the  same  features.  Parameters  as  Fig.  1  but  j=0.14  A. 

AMPLITUDE  BOUNDARY  CONDITION  AVERAGING 

The  modification  of  the  phase  of  the  BC  discussed  above  can  also  be  viewed  as  passing  a  random 
magnetic  flux  through  an  ensemble  of  closed  ring  and  studying  average  properties.  In  this  sense  we  have 
“randomized”  the  locations  of  the  momentum  space  states.  We  have  already  indicated  that  RVA  can 
help  us  transcend  the  limitations  of  this  “bond  phase/magnetic  flux”  approach.  We  can  “randomize” 
electronic  properties  by,  e.g.,  changing  a  local  hopping  or  an  on-site  energy  or  a  Coulomb  repulsion 
somewhere  on  the  chain.  One  could  also  introduce  an  additional  field  and  vary  it  about  zero.  Although 
there  has  been  some  limited  earlier  work  on  using  modified  BCs  in  the  context  of  the  Hubbard  model 


i')),  there  is  as  yet  no  provably  accurate  prescription  for  for  arbitrary  U  and  V.  However,  certain 
intuitive  rules  must  guide  us.  First,  wliatever  change  is  made  to  randomize  must,  of  course,  do  so 
effectively.  Second,  the  change  in  the  system  must  be  negligible  as  the  lattice  size  is  increased  to 
infinity,  for  e.xample,  if  only  one  bond  or  site  is  varied  from  calculation  to  calculation,  then  the  effect 
of  such  a  change  is  immaterial  in  the  thermodynamic  limit.  Finally,  the  behavior  for  small  lattice 
sizes  must  be  illustrative  of  the  infinite-size  limit.  Put  another  way,  one  must  still  be  able  to  make  a 
reasonable  extrapolation  to  the  infinite  chain. 

We  have  found  [2]  an  “amplitude  BCA”  (“scaled- hopping”)  technique  to  be  effective  and  pro¬ 
duce  results  in  good  agreement  with  expectations  based  on  both  strong-  and  weak-coupling  analytic 
arguments.  For  this  method,  we  randomize  by  varying  the  magnitude  rather  than  the  phase  of  the 
“boundary”  hopping  -  I'.e.,  between  sites  1  and  N  -  typicaJly  from  —to  to  -Fto  in  ten  to  twenty  equal 
steps.  The  individual  spectra  are  then  added  together  with  weights  i,  chosen  to  minimize  the  total 
length  of  the  final  curve,  though  giving  each  spectra  equal  weight  yields  virtually  identical  results.  This 
clearly  incorporates  the  special  case  of  using  only  periodic  (x=l)  or  antiperiodic  (i=-l)  rings,  the 
JT/nJT  difference  being  important  for  weak  coupling.  It  also  incorporates  the  case  of  the  open  chain 
(i=0),  which  has  the  “best”  single  BC  size  dependence  (though  with  attendant  “end  effect”  problems). 
Finally,  it  works  in  the  strong-coupling  limit. 

In  Fig.  3c,  we  show  the  spectra  produced  by  the  amplitude  BCA  technique  for  for  a  larger  U  than 
that  in  Fig.  3b  where  phase  BCA  failed.  In  Ref.  |2)  we  showed  that  not  only  are  the  gross  features 
of  the  spectra  obtained  by  amplitude  BCA  in  agreement  with  strong-coupling  calculations  [6],  but  in 
addition  they  show  substantially  more  interesting  detail,  such  as  the  “decoupled-dimer”  peak  located 
at  ~  ff/2  -h  y/{U/2)^  -h  (2to)^  in  systems  with  strong  electron-phonon  and  electron-electron  couplings. 
To  emphasize  the  generality,  in  Fig.  3d  we  show  the  spectrum  obtained  by  varying  U  on  the  first  site 
between  0  and  20  eV.  This  large  U  variation  produces  states  in  the  optical  gap,  but  otherwise  the 
spectrum  is  unchanged  from  Fig.  3c.  We  feel  that  these  results,  and  those  in  Ref.s  [2,3],  demonstrate 
that  RVA  is  an  effective  means  for  reducing  (and  in  some  cases  practically  eliminating)  finite  size 
dependence,  yielding  results  that  can  be  confidently  extrapolated  to  the  infinite  size  limit,  and  which 
are  in  good  agreement  with  known  analytic  results. 
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